The neuron's response at extended timescales 
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Many systems are modulated by unknown slow processes. This hinders analysis in highly non- 
linear systems, such as excitable systems. We show that for such systems, if the input matches 
the sparse "spiky" nature of the output, the spiking input-output relation can be derived. We use 
this relation to reproduce and interpret the irregular and complex 1// response observed in isolated 
neurons stimulated over days. We decompose the neuronal response into contributions from its long 
history of internal noise and its "short" (few minutes) history of inputs, quantifying memory, noise 
and stability. 

PACS numbers: 87.19.11, 87.85. dm, 87.19.1c, 87.19. lr,87.10.Mn,82. 20. Fd 



Many models, especially in biology, are accurate only 
below a certain timescale - due to the existence of addi- 
tional slow processes. If these slow processes are not well 
characterized, it may be hard to predict how they will af- 
fect the dynamics at longer timescales. This is especially 
true if the dynamics are far from equilibrium, highly non- 
linear and contain feedback, a regime where excitability is 
a typical dynamical phenomenon [Ij. Among many types 
of excitable systems (e.g., [Ij and references therein), a 
neuron is a prototypical example - where Action Poten- 
tials (AP - a stereotypical voltage "spike") are generated 
in response to stimulation [2\. AP generation is indeed 
affected by many slow processes [3j - with new processes 
being discovered at an explosive rate [IH^. This may en- 
tail a complex stochastic and history-dependent Input- 
Output (I/O) relation, on multiple timescales [7-9J. In 
general, it is hard to identify, simulate or analyze such 
an I/O due to the large number of processes which are 
unknown or lacking known parameters. 

We find the situation simplifies considerably if we 
use (experimentally relevant [1Q--14J) sparse spike inputs, 
similar to the typical output of the neuron (Fig. [T]A.&B). 
We derive, for a general biophysical stochastic neuron 
model (Eqs. [l]|3| with a few assumptions, a concise de- 
scription for the I/O (Eqs. [6][7|) based on biophysically 
meaningful parameters. This I/O is well described by 
an 'engineering-style' block diagram with feedback (Fig. 
[ip), which can be used to decompose the effects of noise 
and input on the response. Beyond the conceptual lucid- 
ity, such a linear I/O allows the utilization of well known 
statistical tools to derive all second order statistics, con- 
struct linear optimal estimators and perform parameter 
identification. These results hold numerically, even some- 
times when our assumptions break down. 

We demonstrate the utility of our results on recent 
experiments [13j where synaptically isolated individual 
neurons, from rat cortical culture, were stimulated with 
extra-cellular sparse current pulses for a unprecedented 
duration of days. The neurons exhibited 1//" statistics 
[15j, responding in a complex and irregular manner from 
seconds to days. Using our results, we are able to repro- 
duce and analyze the origins of this behavior in a bio- 
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FIG. 1: Mathematical analysis - schematic summary 

A Our aim: to find the I/O relation between inter-stimulus 
intervals (Tm) and Action Potential (AP) occurrences (Ym) - 
for a general biophysical neuron model (Eq. T[3 ). B An AP 
"occurred" if the voltage V crossed a threshold Vth following 
the stimulus. We assume "sparse" stimulations, so Tm ^ 
Tap- C Main result: conversion of a complex biophysical 
neu ron model to a simple linear model with feedback (Eqs. 
6|7 ), where the parameters (F, d, a and w) are linked to the 



biophysical parameters of the full model (Eqs. T][3 ). 



physical model, showing that slow processes span a wide 
range of timescales - with slower processes being "nois- 
ier", due to low ion channel population numbers. The 
model suggests the 1//" statistics of the response orig- 
inates from the long history of internal noise, while in- 
put fluctuations only affect the response on a (relatively) 
short timescale of a few minutes. 

Full Model The voltage dynamics of an isopoten- 
tial neuron are determined by ion channels, protein pores 
which change conformations stochastically with voltage- 
dependent rates [16j. On the population level, such dy- 
namics are generically very well described by models of 
the form [T7H19J 



V = /(F,r,s,/(t)) 
r = A,(F)r + B,(F,r)C, 
s = A,(F)s + B,(F,s)C, 



(1) 
(2) 
(3) 
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with voltage V, stimulation current /(t), rapid vari- 
ables r (e.g., m^n^h in the Ho dgkin- Huxley (HH) model 
[2]), slow variables s (e.g., slow sodium inactivation [20j), 
rate matrices A^/^, white noise processes ^^^^ (with zero 
mean and unit variance), and matrices B^/^ which can be 
written explicitly using the rates and ion channel num- 
bers [21J (D = BB^ is the diffusion matrix [2T, 22j). For 
simplicity, we assumed r and s are not coupled directly, 
but this is non-essential [H]. The parameter space 
can be constrained [19j, since we consider here only ex- 
citable, non-oscillatory neurons which do not fire sponta- 
neously and which have a single resting state - as common 
for cortical cells, e.g., [13j. Such biophysical neuronal 
models (Eqs. T][3) are generally complex non-linear mod- 
els, containing many variables and unknown parameters 
(sometimes ranging in the hundreds [25, 26]), not all of 
which can be identified [27J. Therefore, such models are 
notoriously difficult to tune, highly susceptible to over- 
fitting and computationally expensive [28-30J. Also, the 
high non-linearity usually prevents exact mathematical 
analysis of such models at their full level of complexity 
L31J. 

Model reduction However, much of the complexity 
in such models can be overcome under a well defined and 
experimentally relevant settings [101414]. if we use sparse 
inputs, similar in nature to the spikes commonly pro- 
duced by the neuron. This is done by "averaging out" 
Eqs. [l][3] using similar methods to those in [19j. Specifi- 
cally, suppose / {t) is a pulse train arriving at times {tm} 
(Fig. [iJV, top), so Tm = tm+i-tm > ^Ap with TAP being 
the timescale of an AP (Fig. [l^). Our aim is to describe 
the AP occurrences Ym^ where Ym = 1 if an AP occurred 
immediately after the m-th stimulation, and otherwise 
(Fig. [iJV, bottom). To do so, we need to integrate Eqs. 
[l][3] between tm and tm+i- Since Tm ^ tap the rapid 
AP generation dynamics of (V, r) relax to a steady state 
before tm+i- Therefore, the neuron AP "remembers" any 
history before tm only through = sitm)- Given s^, 
the response of the fast variables (F, r) to the m-th stim- 
ulation spike will determine the probability to generate 
an AP. This probability, pap (s), collapses all the relevant 
information from Eqs. [l][2j and can be found numerically 
from the pulse response of Eqs. [l][2]with s held fixed [24J. 
In order to integrate the remaining Eq. [3] we define the 
averaged rate matrix 

^iXm^Tm) = rApT-i(r^A+ + (l-r^)A_) 
+ (l-TApT-i) Ao, 

where A+, A_ and Aq are the averages of Ag during an 
AP response, a failed AP response and rest, respectively. 
Assuming Tm <C Tg, we obtain, to first order 



Sm+l — -|- Tm-^ iYrri') Tm) H~ I^rn • 



(4) 



T^A {Ym-,Tm))- Note that this simplified linear discrete 
time map has far fewer parameters than the full model, 
since it is written explicitly only using the averaged mi- 
croscopic rates of s (through A and D), population sizes 
(through D), the probability to generate an AP given s, 
Pap (s), and the relevant timescales. This effective model 
exposes the large degeneracy in the parameters of the full 
model and leads to significantly reduced simulation times 
and mathematical tractability. 

Linearization Intrinsic ion channel noise can be ex- 
ploited to linearize the neuronal dynamics, rendering 
it more tractable than the (less realistic) noiseless case 
|19| . Suppose that {Tm} has stationary statistics with 
mean so that tap <C Tm <C Tg with high probabil- 
ity. Since s is slow and AP generation is rather noisy in 
this regime [19j (so pap (sm) is slowly varying [24J), we 
assume a stable excitability fixed point exists, so per- 
turbations = — are small and we can linearize 
Pap (sm) ^ P* + w^s^. The mean AP firing rate can be 
found self consistently (and rather accurately. Fig. [2]/ 
from the location of the fixed point 



{ym) =P* =PAP (s* {p^,T^)) 



(5) 



while the perturbations around the fixed point are de- 
scribed by the linear system 



Sm+i = Fs, 
Y 



dTrr 



oYrr, 



(6) 
(7) 



with F = I + T*A(p*,T*), (n^n^) = T*D (p*, T*, s*), 
em is a white noise process (other parameters in [37J). 
This linear I/O, which contains feedback from the 'out- 
put' Ym to the state variable (Fig. ^p)-, can be very 
helpful mathematically and all its parameters are directly 
related to well motivated biophysical variables. More- 
over, This formulation makes it now possible to construct 
optimal linear estimators for Ym and [32] (Fig. [2^), 
perform parameter identification, and use standard tools 
|22j to find all second order statistics in the system, such 
as correlations or Power Spectral Densities (PSD). For 
example, the PSD for Ym is 

SyU) = wTHe(/)(D(p*,T„s*)+dd^5T(/))Hj(-/)w 



where Wm is a white noise process with zero mean 
and variance T^D (F^, T^, s^) (defined similarly to 



+ T*p4l-p*)|l + w^H,(/)a| 

where He (/) = (27r/z — A (p^,T*) — T~-^aw^) for 
/ <C T~^ . Again, note the large degeneracy here - many 
different parameters will generate the same PSD. Other 
immediately derivable statistics are 5's(/), the cross- 
PSDs SsT if) , SsY if) , Sty (/) and also the respective 
correlations. Our exact results agree very well with 
the numerical solution of Eqs. [l][3] (Fig[2]3-D), even in 
some cases where the underlying assumptions do not hold 
(specifically if Tm ~ Tg and is unstable [24J). 



(8) 
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FIG. 2: Comparing the mathematical results with the 
numerical simulation of Eqs. [TKsl for the stochastic 
HH model with slow sodium inactivation [191 [24| with 
Tm = n = 50ms. A Firing probability (T"^) (Eq. Isl) 
for different currents (/q = 7.5,7.7,7.9,8.1,8.3 /xA from bot- 
tom to top). B Optimal linear estimation of s. C-D Power 
spectral densities Sy (/) and Ss (/). 'Map' is a (10^ faster) 
simulation of Eq. [4] together with pap (sm), while Approx' 
is the analytical expressions (e.g., Sy is given by Eq. |8|. 
lo = 7.9/iA in B-D. For more complex models see [24J. 



Modeling the Experiment Next, we use the ex- 
pression for the PSD (Eq. |8| to find conditions un- 
der which the experimental results of [13j can be repro- 
duced, given our assumptions. We will show that kinetic 
processes must possess a large range of slow timescales, 
where the number of relevant ion channels with each 
timescale scales with an exponent a. Fitting a specific 
model we reproduce the experimental results in Fig. |3] 

Previous work [T^ used a stochastic HH model with 
slow sodium inactivation to reproduce the experimental 
results of [13] up to a timescale of minutes. However, 
this model cannot explain dynamics on longer timescales 
(Fig. [7p). Specifically, we require Sy (/) ^ for / < 
IQ-^Hz, with a ^ 1.4 when = T*, as in [13]. If the 
total number of ion channel states of all channel types is 
finite, then we can decompose Eq. [8] 



Sy if) = Co 



M 

E 

i=l 



CiXi 



where q are some (derivable) constants and the poles 
are solutions of the characteristic equation 



|AI-A(p*,T*)-T, 



= 0. 



In order to approximate a PSD of the form /~", we re- 
quire the poles to cover a large range |15j. Though the 
T~^aw^ feedback term can tune the location of the poles 



(through the variable a, see Fig. [4p), comparison with 
experiment implies that the observed dependence 
was not generated in this way, since it exists even near 
the critical stimulation frequency, where 1, w ^0. 

Therefore, the eigenvalues of A (p*,T*), the average rate 
matrix, must span a large range of (inverse) timescales. 
However, the existence of this range is not sufficient - we 
also require [15] q oc A^^" so that 
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Analyzing Eq. [8] we find that we can generate this in a 
robust way (independent of T^), that is consistent with 
other experimental observations, by having a scaling rela- 
tion in Ni, the number of the corresponding ion channels 
(affecting D(p^,T^,s^), the variance of in Eq. [6|. 
Specifically, by setting Ni cx A~", implying that slower 
processes are noisier. To fit a specific model consider 
the stochastic HHS model which augments the basic HH 
model with a slow sodium inactivation variable |T9j [24] . 
We extend this model by replacing the sodium inactiva- 
tion variable s with ^i/^ •> where si is identical to 
5, and for {si}^^ rates scale as and the channel 
numbers scale as Ni cx e"*. In order to obtain a = 1.4 
for the duration of the experiment we require e = 0.2 
and M = 5. This reproduces well the observed scaling 
relations (Fig. [3|. 

Predictions First, we consider the duration of the 
neuronal memory. To quantify this more precisely, we 
note that in the frequency domain we can use spectral 
factorization [32j and write the response of the linear 
system (Eqs. [6][7|) as 

Y (/) = (/) f (/) + ifnoise (/) V if) (9) 

with i^signai (/) = w^Hg (/) d and u (/) is the Fourier 
transform of a zero mean white noise representing in- 
ternal neural noise. In our model, the behavior is 
generated by fluctuations in internal neuronal noise, so 
Hnoise if) ^ However, i^signai (/) c for / ^ 

since it is not affected by A^^. This entails that the neu- 
ron possesses a long memory for its intrinsic fluctuations 
but a "finite memory" of its input - i.e. perturbations in 
Ym due to perturbations in Tm will decay exponentially 
with a finite timescale. Specifically, in the fitted model, 
^signal (/) has the shape of a simple low pass filter with 
a timescale of ~ 100 sec. This could also be probed ex- 
perimentally directly by applying a sinusoidal Tm input 



Tm = T,^^Ti sin (2^/,T*m) , 



(10) 



since the linear response of the model will generate a 
direct probe of i^signai (/) (Fig. |4|V). Note that this pre- 
diction of input memory timescale of ~ 100 sec is valid 
only for spiking input, in the context of our model. In 
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FIG. 4: Input memory in fitted model. (A) 

'Approx':|i7signai (/)| of the fitted model has a low-pass fil- 
ter shape with cutoff at / = 10~^ Hz. We predict this shape 
could be probed directly by peaks in Fourier transform of Ym 
('simulation'), if sinusoidal input (Eq. 10) is used. (B) Same 
model when stimulated with 1// stimulation pattern in com- 
parison to experiment in frequency ranges where high SNR 
allows reliable estimation of Sty (/)• 



FIG. 3: The measures of "scale free" rate dynamics in 
the extended stochastic HHS model - comparison of the 
experimental data from [13 and a simulation of the extended 
HHS model (solid and dashed lines, respectively). We use 
here the same measures as in Fig. 6 in [13 : A The firing 
rate fluctuations estimated using bins of different sizes (10, 
30, 100, 300 sec) and plotted on normalized time axis (units 
in number of bins), after subtracting the mean of each series. 
B CV of the bin counts, as a function of bin size, plotted on 
a log-linear axis. C Firing rate periodogram. D Detrended 
fluctuations analysis. E Fano factor (FF) curve. F Allan 
factor (AF) curve. G Length distribution of spike-response 
sequences, on a half-logarithmic axes. H Length distribution 
of no-spike-response sequences, on a double-logarithmic axes. 
For additional details see and references therein. 



fact, continuous input signals may very well generate long 
memory effects [33j. 

Second, we quantify the "noisiness" of the neuronal re- 
sponse. In the fitted model the neuronal response under 
periodic stimulation is very noisy - in the sense that linear 
optimal estimation of performs similarly to the trivial 
predictor (Ym = \{Ym) — 0.5]), even if is fully known 
(since has a large variance). In order to improve pre- 
dictability of Ym, we can increase the variability of the 
input Tm- To test this we examined data from a similar 
experiment where the variability of Tm was higher than 
the internal noise at certain frequencies, so it was possible 
to estimate Syt if) = ^signal (/) St (/) in these regions 
(which is generally hard if internal noise is high [34J), 
which seems to correspond well with our fitted model 
(Fig. [4^). Therefore, input variability in inter-spike in- 
tervals Tm can increase signal to noise ratio, as was pre- 
viously observed for current amplitude variability over 
shorter time scales of milliseconds [35j. 

Finally, consider the stability of the neuronal response 
for even longer timescales. Generally, Sy (/) ~ as 
/ ^ with a > 1 implies Var (Ym) ~ m"~^ as m ^ 00. 
However, Var (Ym) < 0.25 since Ym is binary, and there- 



fore the scaling behavior must have some cutoff. Extrap- 
olating the experimental results suggests that for a = lA 
this cutoff could be reached approximately only after four 
years of ongoing experiments - which is comparable with 
the lifespan of a rat (in comparison, a = 2 gives a cutoff 
of a few days). In general, M = 5 poles should allow 
coverage of all this timescale range [15j. However, in our 
model, the scaling is limited by channel numbers. Since 
typically 1 < A^i < 10^ [36j, scaling in Ni can generate 
in the PSD over 6/ a Orders of Magnitude (OM). 
For a = 1.4, we get 4.2 OM, while [13j covered about 3 
OM. 

To summarize this section, we used the linearity of 
our derived I/O to decompose the contributions of inputs 
and internal noise to the response of the neuron. This 
decomposition shows that even though the neuron can 
"remember" its intrinsic fluctuations over timescales of 
days, its memory of past inputs can be relatively "short". 
For example, the input memory of our fitted model decay 
exponentially with a timescale of ~ 100 sec. We suggest 
an experiment to test this directly. Additionally, this 
linear decomposition allows us to quantify the signal-to- 
noise ratio of the I/O and show that it increases with 
input spike time variability Finally, we set upper limits 
on the timescale range of the observed 1//" behavior. 
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